2, 2008 5:54 WSPC/INSTRUCTION FILE iHMC-ijmpc-rev 



International Journal of Modern Physics C 
© World Scientific Publishing Company 



PROPERTIES OF ITERATIVE MONTE CARLO SINGLE 
HISTOGRAM REWEIGHTING 



MARTIN GMITRA 

Department of Theoretical Physics and Astrophysics, P.J.Safdrik University, 
Park Angelinum 9, 040 01 Kosice, Slovak Republic 
gmitra@kosice.upjs.sk 

DENIS HORVATH 

Department of Theoretical Physics and Astrophysics, P.J.Safdrik University, 
Park Angelinum 9, 040 01 Kosice, Slovak Republic 
horvath. denis&gmail. com 

Received Day Month Year 
Revised Day Month Year 

We present iterative Monte Carlo algorithm for which the temperature variable is at- 
tracted by a critical point. The algorithm combines techniques of single histogram 
reweighting and linear filtering. The 2d Ising model of ferromagnet is studied numer- 
ically as an illustration. In that case, the iterations uncovered stationary regime with 
invariant probability distribution function of temperature which is peaked nearly the 
pseudocritical temperature of specific heat. The sequence of generated temperatures is 
analyzed in terms of stochastic autoregressive model. The error of histogram reweight- 
ing can be better understood within the suggested model. The presented model yields 
a simple relation, connecting variance of pseudocritical temperature and parameter of 
linear filtering. 
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1. Introduction 

The critical phenomena, typical of many body systems, have attracted remarkable 
scientific interest since a long time. The Monte Carlo (MC) stochastic techniques^ 
serve to estimate thermodynamic averages from samples of the configuration space. 
Many advanced simulation approaches are rooted in elementary updates of Metropo- 
lis algorithm El which generates a Markovian chain of spin configurations. 

Most of the efforts try to improve MC method within the critical region. The 
efficiency of MC method have been increased by the development of sampling and 
data processing te chni ques. The reweighting technique applied anew to classical 
statistical systems enables to perform a random sampling for one distribution 
function, and to calculate quantities of interest for similar distributions. To be more 
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concrete, let us assume that a sample of N configurations (g = 1,2,..., N) of quan- 
tity Q have been accumulated for a fixed temperature Tj. The samples Qq are used 
to construct a continuous temperature T dependence of a single histogram approx- 
imation {Q)T.Tt of the canonical thermodynamic average {Q)t- The interpolation 
called " reweighting on the fly" E! ig based on the formula 



where Eq is the energy (in Boltzmann constant units) of configuration q. Because of 
the finite length of the MC run, Eq. ^ provides reliable results only for a relatively 
narrow range of T values around Tt . When T differs too much from Tt , it causes 
an increase of the statistical errors and may lead to unreliable results. The reliable 



in Ref. 8 has shown that reweighting causes a systematic shift in the height and the 
location of the specific heat peak, which depends on the relative position of Tt and 
T. It turns out that the shift decreases as ~ iV^^/^ only for N much larger than the 
number of degrees of freedom. Systematic errors, eventually affected by sampling 
autocorrelations, may also take place at finite N. The first quantitative study of the 
statistical errors in reweighted dataElhas demonstrated that theoretical calculation 
of the error is both difhcult and time consuming so it may not be justified for all 
studies. In this paper we present rather general iterative treatment based on the 
single histogram reweighting which implicitly takes into account both the systematic 
errors of reweighting and the statistical errors to generate stochastic fluctuations 
near a critical point. 

A well known systematic error, present in the histogram method, is an overes- 
timation (underestimation) of {E)T.Tt when T -^Tt {T ^ Tt), respectively 1^. This 
fact and the general idea of feedback regulation!^ are used to construct an iterative 
algorithm. The aim of the suggested method is to estimate the critical temperature 
or other parameters of interest, as well as to perform an analysis and partial elimi- 
nation of their errors. The method combines single reweighting, linear filtering and 
an iterative treatment. The incorporation of iterations is inspired by a general class 
of probabilistic models The iterations generate data that are further analyzed 
in terms of a suggested autoregressive moving average (ARMA) model. 

The single-spin flip Metropolis algorithm is applied here in order to obtain av- 
erages for moderate lattice sizes. The cluster methods are not used due to the 
efficiency reasons discussed in Ref. El The efficiency of the algorithm should also 
be increased by hybrid schemes 'i^. However, because of many side effects that can 
show up, a separate study of hybrid schemes is needed. 

The paper is organized as follows. In Sec.|21the iterative self-organized algorithm 
is suggested. In Sec. 01 the autoregressive phenomenological model of the stochastic 



{Q)t,t, 




(1) 



range of T values also decreases as the system size increases 



'67 



. A qualitative study 



2, 2008 5:54 WSPC/INSTRUCTION FILE iHMC-ijmpc-rev 



Properties of iterative Monte Carlo single histogram reweighting 3 

iterative search process is proposed and the corresponding relations for second- 
order statistics are derived. The model is used to parametrize MC data obtained 
for ferromagnetic 2d Ising spin models on L x L square lattices (see SecQJ. 

2. The iterations near the criticality 

We start by defining a process where temperature is driven by some response func- 
tion. The stochastic process can be formally described by the second-order recursive 
formula 

T^I\=Ft,.n, (2) 
Tt+i = vTli\ + (1 - v)Tt . (3) 

It generates sequence { Tt,Tt=o — To,t — 1,2, . . . , iVovcnt}- The stochastic function 
F : Tt ^ T^^\, which provides a preliminary histogrammatic estimate Tj^i of 
unknown critical temperature T* , is defined implicitly by 

F : C'rpMs Tt N = niax C'T.Tt,N ■ (4) 

*■ Te(0,oo) 

The key point of the reweighting is to localize the extremal value of the response 
function C evaluated for N MC samples. Commonly, the response function can be 
defined by the first and second moments of energy or magnetization, respectively. 
The role of filtering, defined by the plasticity parameter 77 in Eq. is to weaken 
fluctuations ^1 of the T^^\ sequence. 

Let us assume the invariant probability distribution function (pdf) of the form 

I 

p(T) = 5T,Tt = lim ^ Sr^Tt , (5) 

where S is the Kronecker's symbol and is the arithmetic average (which differs 
from standard thermodynamic average (. . .)). From now on T* is associated with 
pseudocritical temperature Tc{L) of the specific heat. For Gaussian-like p{T) we 
suppose 

T*~ lim J2Tt^Tt (6) 



/ dTp{T)T. 
Jo 



As mentioned above, the specific heat can be estimated within a histogram reweight- 
ing approach using the fluctuation-dissipation relation 

At the thermodynamic limit N ~^ 00, the iterative process becomes deterministic 
and tends to an unique fixed point T* — Ft* ^n^oc- For finite N, a convergence 
criterion have to be formulated in a pure statistical sense. The variance 



Vt = {Tt-T*f (8) 
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has been chosen as a suitable measure that reflects the convergence. 
3. The model of the iteration process 

The iterative process, described by Eqs. lO, |@J supplemented by Eq. lO, 
consists of nonlinear and stochastic rules. As an auxiliary tool of analysis, let us in- 
troduce a phenomenological model to parameterize the simulated statistics. Within 
the model, the explicit form of the stochastic function [see Eq. is 

FT,^N = aT* + {l-a)Tt+^t , (9) 

where an additive noise S,t describes uncertainty. The convergence rate towards T* 
is determined by the parameter a, which contains information of the systematic 
reweighting error. From Eqs. Q and Eq. ^ follows that 

Tt+i ^ arjT* + {1 - a7])Tt + v^t , (10) 

where the noise term is controlled by the coefficient 77. The model discussed here can 
be considered as an application of ARMA(1,1) model The solution of Eq. (|10|) 
can be written simply as 

t-i 

Tt^T* + (To - T*)B-' + i]B J2 ^qB"-* , (11) 

q=Q 

where B = 1/(1 — arj). The above formula, written for any dependence, is general. 
At this point, let us consider the statistical properties of Gaussian noise for 
specified by 

cr = o, (12) 



66' =ASt,t' , (13) 

where A > 0. Within the limit of a second order statistics in Tt, from Eq. (|ll|l and 
Eq. follows 

Tt ^ T* + [Tq - T*)B-^ . (14) 

Therefore, the stationary regime is defined by the invariance Tt — T* , which stays in 
agreement with Eq. . If the solution given by Eq. 1)11(1 is substituted into Eq. , 
one obtains 

Vt = (To - T*f + 1^ ( 1 - ) , (15) 

which is equivalent to the recursion Vt+i = Arj^ + Vt B^^. Terms proportional to 
B^^* describe a transient regime of the initial temperature Tq, which the system 
forgets exponentially fast. Terms [Tq ~ T*)^^-^* and Ar]'^B'^''^-*^ /{I - B'^) are de- 
terministic and stochastic contribution to Vt , respectively. The variance Vt remains 
finite for t — > 00 if > 1, that is if < ayy < 2. The transient regime can be 
characterized by 

Vi= Ajf + {l-arif{To-T*)\ (16) 
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whereas, to obtain the stationary regime, one needs to calculate 

It results that Voo diverges at a ^ 2/77 or equivalently B ^ —1. The characteristic 
transient time, 

ru- = -^^ 1' (18) 

In |1 — ar]\ 

follows from its definition — exp(— 2t/rtr). The transient time has a singu- 

larity for 77 (l/a)^ and B — > ±00. In this limit the model describes immedi- 
ate convergence (rtr 0) and variance Voo — Ajo?. A singular parametric line, 
1 — arj = 0, splits region B^ > 1 into alternating {B < —1) and unidirectional 
convergence (B > 1) of Tt near T* . 

Regarding the convergence properties, it is useful to introduce 77 = 77^ as the 
value of the highest initial convergence rate, defined by the condition dVi/d77| ^=,,^ = 
0. This yields 

a 1 



For this extremal value is Vi(77ni) — {A/a)ritn- 

Finally, we introduce some supplementary remarks about the limit of the van- 
ishing error. Clearly, for iV — > 00 we assume a ^ 1~ . The consequence of Eqs. (II 7|) 
and lfTH|l is 77 — > 0+, that implies Vt^oo 0. Unfortunately, because of rtr 00 
this cannot be attained in practice. 



4. Numerical results 

In order to explore the proposed model by numerical simulation, we chose the 
2d Ising model with the Metropolis algorithm as testing ground. Most of the data are 
accumulated on square lattices of linear dimension i = 10 and periodic boundary 
conditions. The methodology is expected to be quite general and of straightforward 
application for both discrete and continuous systems. 

Every iteration step t = 1, 2, . . . , A'^cvcnt of algorithm, based on Eq. and 
Eq. consists of the following main blocks (i)-(iv): 

(i) equilibration that includes 10^ MCSS (Monte Carlo steps per spin) for constant 
Tt. 

(ii) accumulation of N pairs {Qq,Eq) into histograms. In order to reduce autocor- 
relations, successive samples are separated by one MCSS. 

(iii) localization of the maximum Crphis^ jy by a gradient algorithm. Recall, that C 
is computed from Eq. (P) and Eq. (0) for (E)T,Tt and {E'^)T.Tf 



2, 2008 5:54 WSPC/INSTRUCTION FILE iHMC-ijmpc-rev 



6 M. Gmitra, D. Horvdth 

(iv) prediction of next Tt+i using the filtering given by Eq. 

Fig. ^ illustrates the properties of transient regimes, which are affected by the 
choice of Tq and rj. The T^-dependences of Vi are presented in Fig. |2a). One may 
see that the agreement of the simulated data with Eq. (|16|l is quite satisfactory. 
Typical values of the parameters determined by the fit are gathered in Tabled 




To/T,{L) 

Fig. 1. The transient stages of the iterative regime characterized by the minimum number of 
iteration steps n needed to satisfy the stop criterion \ T„ — Tc{L) \ < 0.01 for resulting T* = 
Tc{L) = 0.586141. The number n is plotted as a function of initial temperature Tq [distinct from 
Tc(L)] for several r/. Calculated for fixed = 10^. 



Table 1. The transient regime param- 
eters for To = 0.6 obtained by fitting of 
Eq. JTO}. 



N 


a 


A 






rim 


10^ 


1.01805 


1.398 ■ 


10- 


4 


0.5770 


10* 


1.00281 


1.227 ■ 


10- 


5 


0.9376 


10^ 


0.99937 


1.282 ■ 


10- 


6 


0.9939 


10"^ 


0.99974 


1.216 ■ 


10- 


7 


0.9996 



For N oo, we observe that rj^i — > 1^, which is in agreement with our ex- 
pectations. The insert of Fig. [^b) illustrates the dependence of ?7m(7o) with the 
special cases Tq = Tc{L) and Tq ^ Tc{L), which clearly coincide with Eq. H19|l . 
Additional numerical study reveals the discrepancy for Tq <C Tc{L), that is caused 
by a non-trivial decrease of a by decreasing Tq {a < 1). The A^-dependence de- 
picted in Fig.[2Ib) is constructed for an initial value Tq = 0.6. There is a certain N 
value for which a = 1 and iterations cross from alternating to unidirectional conver- 
gence. In other words, statistical errors compensate systematic ones. In agreement 
with Ref. IHl we have verified statistical errors of MC averages to be proportional 
to 1/VN, which implies A^l/N. Actually, the fit provides A~0.13/N. 
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0.001 




ICp 1^ up 10^ 

N 

Fig. 2. The analysis of transient regime: (a) The rj-dependence of variance Vi obtained for initial 
To = 0.6. The lines correspond to the fits of r;-dependences via Eq. I16i . The parameters are 
collected in Tabled] (b) The a{N) dependence for To = 0.6. The inset depicts r)m(To) for A? = lO". 

If the number of iterations is much larger than rtr, the initial conditions are 
forgotten and thus invariant pdf 's of temperature variable are formed (see inset in 
Fig.EJl. The calculation of T* = Tc{L) {t > n,. > 1) is controlled by 77. We have 
verified VV* oc rj/^/N for the stationary regime. The error increases monotonically 
with 77 for a given N as in Fig. |21 We have compared the simulation and analytical 
77-dependences of Voo for A and a, given in Tabled and those determined by a fit 
of Eq. H17|) for N — 10^. The comparison reveals a discrepancy of the A parameter 
determination. The A parameter, determined by the fit of Eq. (|17|l to numerical 
data, increases its value about 6.2 x 10^*. A further extension of the model, e.g. 
taking into account a weak colored noise, should eliminate this distinction. 

According to Eq. JHl, the mean values of invariant pdf's (see inset in Fig. 
were used to estimate Tc(lO) = 0.586141. The method is straightforward applied 
for L = 10, 20, . . . , 100. Subsequently, the values Tc{L) are interpolated by the 
finite size formula Tc{L) = Tc + h/L. This yields to estimate the thermal coeffi- 
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cient b = 0.1889 ± 0.0033 and the critical temperature as = 0.5673 ± 0.0001. 
The last estimation stays in a satisfactory agreement with exact value T^^ = 
[2 In (l + V2)]^^ = 0.56729. The problem of critical indices can be handled ef- 
fectively in the frame of a two-lattice iterative algorithm, too^l In that case, the 
objective function given by Eq. Q has to be written in terms of the Binder cumu- 
lant, F : C/ti- ,t,,7V 



min \UTTt n{L2) — Utt, n{Li)\, for two different system 

Te(0,oo) ' ' ' ' 

sizes. According to the fact, that histogram reweighting allows a superior determi- 
nation of the derivation U!p^^_^ Ni-^) = ^UT,Tt,N{L)/dT\T=Tt+n we considered, as 
an example, the exponent v of the correlation length 



In 



^Tt+i,Tt,jv(-^l) 



Within the proposed algorithm, we generated a sequence of estimates (for each 
iteration step separately) of the critical exponent for the lattice pair 2Li = L2 = 20. 
Considering the peak of nonsymmctric pdf of {l/i')t we determined l/v''^ = 1 with 
a relative error of 1.35%. 




Fig. 3. The 77-dependence of Voo compared with Eq. 1171 for parameters A, a and N = 10^ 
given in Table [Tl The error bars are calculated using 5000 independent runs. The inset depicts a 
stationary pdf of Tt constructed for Nevent = 10^ and 77 = 1. 



5. Conclusions 

An iterative process that is useful around the critical point and includes Metropolis 
single spin-flip algorithm, single- histogram reweighting and linear filtering is studied 
numerically. The properties of this iterative algorithm are studied for the 2d Ising 
ferromagnetic model on a square lattice. The total errors in the estimate of pseu- 
docritical temperature are obtained. We demonstrated how the accuracy can be 
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affected by rj. However, one should be careful with an extremely small rj to avoid 
an enormous increase of transient time. Our analysis demonstrates that it is insuffi- 
cient to use several reweighting trials to estimate total errors. The most interesting 
iterations are of this kind that compensate the systematic errors introduced by the 
histogram reweighting approach. Despite of the success of the linearized model, 
which allows interpretation of numerical data, further improvement of the proposed 
methodology is still expected. For instance, one should take into account non-linear 
corrections in Eq. (jlOfl and/or the colored noise (that is, with autocorrelations of 
error terms) in Eq. (|13|l . We believe that nontrivial dependences of a and 77^ (see 
e.g. inset of Fig.[2tb)) could be further explained considering microscopic expression 
for the density of states'^. 

Acknowledgments 

The authors would like to thank the Slovak Grant agency VEGA (grant 
No. 1/2009/05) and the agency APVT-51-052702 and the internal grant of the Fac- 
ulty of Sciences of Safarik University VVGS 2003 for financial support. The authors 
express their thanks to unknown referee for valuable corrections and suggestions. 

References 

1. K. Binder, D.W. Heermann, Monte Carlo Simulation in Statistical Physics (Springer, 
Berlin, 1998). 

2. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, J. Chem. Phys. 21, 
1087 (1953). 

3. Z.W. Salsburg, J.D. Jacobson, W. Pickett, W.W. Wood, J. Chem. Phys. 30, 65 (1959); 
A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988). 

4. A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989). 

5. A.M. Ferrenberg, D.P. Landau, R.H. Swendsen, Phys. Rev. E 51, 5092 (1995). 

6. A.M. Ferrenberg, D.P. Landau, Phys. Rev. B 44, 5081 (1991). 

7. D.P. Landau, Journ. Magn. Magn. Mater. 200, 231 (1999). 

8. E.P. Miinger, M.A. Novotny, Phys. Rev. B 43, 5773 (1991). 

9. L.P. Kadanoff, Physics Today (March 1991), p. 9. 

10. H. Robbins, S. Munroe, Ann. Math. Stat. 22, 400 (1951). 

11. N. Ito, G.A. Kohring, Int. J. Mod. Phys. C 5, 1 (1994). 

12. J.A. Plascak, A.M. Ferenberg, D.P. Landau, Phys. Rev. EQ5, 066702 (2002). 

13. J.C. Principe, N.R. Euliano, W.C. Lefebvre, Neural and adaptive systems: Fundamen- 
tals through simulations (John Wiley & Sons, Inc. 2000). 

14. P.J. Brockwell, R.A. Davis, Introduction to Time Series and Forecasting (Springer, 
New York, 1996). 

15. D. Horvath, M. Gmitra, Z. Kuscsik, Czech. J. Phys. 54, 921 (2004); 
D. Horvath, M. Gmitra, Int. J. Mod. Phys. C 15, 1269 (2004). 

16. P.D. Beale, Phys. Rev. Lett. 76, 78 (1997). 



